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Abstract. This text proposes a fast, rapidly convergent Nystrom method for the solution of the 
Lippmann-Schwinger integral equation that mathematically models the scattering of time-harmonic 
acoustic waves by inhomogeneous obstacles, while allowing the material properties to jump across 
the interface. The method works with overlapping coordinate charts as a description of the given 
scatterer. In particular, it employs “partitions of unity” to simplify the implementation of high- 
order quadratures along with suitable changes of parametric variables to analytically resolve the 
singularities present in the integral operator to achieve desired accuracies in approximations. To deal 
with the discontinuous material interface in a high-order manner, a specialized quadrature is used 
in the boundary region. The approach further utilizes an EFT based strategy that uses equivalent 
source approximations to accelerate the evaluation of large number of interactions that arise in 
the approximation of the volumetric integral operator and thus achieves a reduced computational 
complexity of 0{N log N) for an Y-point discretization. A detailed discussion on the solution 
methodology along with a variety of numerical experiments to exemplify its performance in terms 
of both speed and accuracy are presented in this paper. 


1. Introduction 

The solution of direct acoustic scattering problem where the goal is to obtain the scattered wave 
produced as a result of an interaction between a given incident wave and a given bounded inho¬ 
mogeneity continues to constitute one of the most challenging problems in computational science, 
especially where applications of interest require such computations to be carried out for large pen¬ 
etrable scatterers with complex geometries. Formally, the d-dimensional direct acoustic scattering 
problem that we consider in this paper for d = 2 and 3 is described as follows: given an obstacle fl, 
a bounded open subset of M'^, with a smooth boundary dQ, and an incident time-harmonic acoustic 
wave u* satisfying 

(1) Au*(a;) -I- K^u\x) = 0, x e 

where k = uj/cq is the wavenumber, u is the angular frequency, and cq is the constant speed of 
wave outside the inhomogeneity 0, find the total acoustic field u that satisfies [T] 

(2) Au(®) -b n^n^{x)u{x) = 0, ® G 

with the refractive index n{x) = co/c{x), where c, the speed of acoustic wave, is allowed to vary 
with position within Q and the scattered field := u — u* satisfies Sommerfeld radiation condition 

(3) lim f - iku^^ = 0, 

r^oo \ or J 

where r = \\x \\2 = 


1 
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One of the main challenges in obtaining numerical approximation of the solution to the scattering 
problem arise from the need to accurately describe highly oscillatory functions in large domains. 
As a certain fixed number of points is necessary to resolve a wavelength, we typically require a large 
computational grid to obtain any meaningful solution and any effort to further reduce the compu¬ 
tational errors demands that the size of discretization grows proportionally in all dimensions. It is, 
therefore, highly desirable to devise numerical methodologies that are efficient as well as high-order 
accurate. In this direction, several methodologies, based on diverse formulations of the problem, 
ranging from differential equation to variational form and to integral equation, have been proposed 
by various authors. In the solution strategies that are based on differential equation or their cor¬ 
responding weak formulations laEllllIlEllIlElElllo], a relatively large computational domain 
containing the scatterer must be used, together with appropriate absorbing boundary conditions on 
the boundary of the computational domain. Thus, these procedures often require a large number of 
unknowns and, hence, lead to large linear systems. In addition, accurate absorbing boundary con¬ 
ditions with efficient numerical implementations are quite difficult to construct; the error associated 
with such boundary conditions typically dominates the error in the computed solution. 

In contrast, the integral equation approach, where the mathematical formulation directly ensures 
that the solution satisfies condition Q by suitably employing the radiating fundamental solution, is 
free from considerations mentioned above, and consequently, does not require solution strategies to 
discretize outside the inhomogeneity. We, therefore, base our numerical treatment of the scattering 
problem on an equivalent integral equation formulation which is given by the Lippmann-Schwinger 
equation mm, 


( 4 ) 


u{x) + K^ / Gn{x,y)m{y)u{y)dy = u\x), x£W^, 


where 

( 5 ) 


in 


’ I exp(iK|a; — ^|)/47r|a; — y|, in M^, 


is the radiating fundamental solution of Helmholtz equation in the free space and m{x) = l — n^{x). 

In recent years, a lot of progress has been made toward numerical solution of the Lippmann- 
Schwinger equation; for example, see [iziiisiiiaiiaiiiiiiziiisiiigiEoiEiiisaEalaEaEeiET]. 
Most fast numerical schemes among these, though high order accurate for smooth scattering media, 
exhibit only linear convergence in the presence of material discontinuity [laiiaiiTlEaEelEZ]. For 
example, in [26], Duan and Rokhlin, introduced an efficient high order quadrature formula which 
utilizes a corrected Trapezoidal rule in conjunction with FFT for fast and accurate approximation 
of volume integration in Q. The high order convergence of the scheme, however, requires the 
refractive index n{x) to be globally smooth. More recently, a high order direct solver with 
computational complexity has been proposed in m- Again, while this algorithm is robust and 
provide accurate approximation for smooth scatterer, it does not exhibit rapid convergence in the 
presence of a discontinuous material interface. A couple of fast techniques, one that rely on the 
use of “discontinuous FFT” [24| for efficient computations while the other uses FFT for accelerated 
evaluation of convolution in polar coordinates after a suitable decomposition of the Green’s function 
via addition theorem m, however, achieve second order accuracy for discontinuous scattering 
configurations. Another approach for solution of acoustic volumetric scattering problem introduced 
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in mm, though high-order convergent, is designed to be computationally efficient only for “thin” 
scattering configurations. This scheme gains high-order convergence through a combination of 
changes of parametric variables (in order to resolve the singularities present in the Green function) 
and by suitably employing “partitions of unity” to yield smooth and periodic integrand away 
from vicinity of target points. Our present approach, in fact, is a non trivial extension of ideas 
presented in mm wherein we obtain a solver for general scattering configurations that exhibits 
computational complexity of 0{N log N) with respect to the grid size N while retaining high-order 
accuracy even in the presence of material discontinuity. We believe that this algorithm is the first 
fast integral equation solver which provide high-order convergence for full volumetric scattering 
problem with a discontinuous material interface. 

The rest of the paper is organized as follows. In Section we present main algorithmic com¬ 
ponents of our numerical scheme. We then present, in Section a detailed account of our ap¬ 
proximation strategy for the integral operator arising in the context of two dimensional volumetric 
scattering problem. The last part of this section presents a series of numerical experiments to 
exemplify the performance of the proposed method, both in terms of accuracy and in terms of 
computational efficiency. This is followed, in Section]^ by a brief discussion on numerical solution 
of corresponding three dimensional scattering problems using a straightforward extension of our 
two dimensional solver and demonstrate its effectiveness via some numerical experiments. Finally, 
our conclusions are summarized in Section [5j 

2. Principal components of the method 

As mentioned earlier, we base our numerical strategy on solving the linear system arising out 
of a Nystrom discretization of the Lippmann-Schwinger integral equation. Given the denseness 
of resulting linear systems where use of direct linear solvers can be prohibitively expensive, we 
employ the matrix-free version of the iterative solver GMRES [28], in fully complex arithmetic, 
for the solution of the discrete form of Q. In this context, we focus our presentation on the 
computational technique for accurate and efficient evaluation of the integral operator 

(6) JC[u\{x) = j G^{x,y)m{y)u{y)dy. 

Q 

We start by describing the scatterer 11 through a collection V = {Vk}k=i ^ overlapping 
coordinate patches where the fe-th patch is homeomorphic to an open set Tik C (0,1)'^ via a smooth 
invertible parametrization ^k- This, in conjunction with a partitions of unity (POU) subordinated 
to the covering V, that is, functions cok{x) : k = 1,..., K, satisfying 

K 

''^^ujk{x) = 1 for all X £ id, 
k=l 

where for each k, iOk £ C°°{Q) with its support contained in Vk, reduces the evaluation of (§ 
to computation of integrals over these K patches. Note here that for all those patches Vk whose 
closure does not intersect with the boundary of il, the corresponding Uk vanishes to high order 
along with all of its derivatives on the boundary of the patch. On remaining patches, however, 
where one of the edges coincide with the boundary of il, ujk clearly does not vanish and, in fact, 
attains the value 1 (see Figure for examples). For the clarity of presentation of the proposed 
quadrature, and to distinguish between patches based on this criterion, we introduce two index 
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Figure 1. Partitions of unity. 


sets, namely, X/ = {fe | n dfi = 0} and = {/c | n dQ / 0} corresponding to interior 
and boundary patches respectively. Now, clearly, the integral can be reexpressed as a sum of 
integrals over boundary and interior patches 

(7) /C[ul(rt) = / G^{x,y)m{y)u{y)ujk{y)dy +/ Gi^{x,y)m{y)u{y)u}k{y)dy. 

keig k&Xi 

In view of this, we introduce a set of quadrature points on the k-th patch, say ‘Ik, by taking the 
image of the equi-spaced grid {{h/Ni,.. ■ ,id-i/Nd-i,id/Nd) \ 0 < ii < Ni,i = 1,... ,d} in [0, 1]'^ 
under the map Further, the union of these grid points, that is, 

K 

(8) X = U ‘Ik 

k=l 

where 

(9) Ik = {^k{ii/Ni, ■ ■ ■, id-i/Nd-ijid/Nd) | 0 < < Ni, £ = 1,... ,d} 

dehnes the complete set of Nystrom discretization points where we seek to approximate the solution 
of equation Q. 

Indeed, the performance of our iterative solver hinges on our ability to approximate the integrals 
in Q accurately, in a computationally efficient manner. In the next section, we provide a detailed 
description of a fast, high order integration scheme in two dimensions that can be readily extended 
to three dimensions, which we briefly outline in Section]^ 

Remark 1. In this text, we often refer to x in equation ^ as a target point whereas points y 
therein have sometimes been called source points. 

Remark 2. While discussing the approximations on boundary patches, we identify the d-th co¬ 
ordinate variable in [0,1]“^ with the transverse parameter. We, therefore, use the notation t = 
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(ii,..., td-i,td) to denote a point in the corresponding parameter space and assume that the bound¬ 
ary coincides with tj, = 0. 


3. Fast and accurate evaluation of integrals in two dimensions 

The difficulty in accurately computing integrals in Q is significantly more when the target point 
X lies in the integration patch Vk compared to the case when it does not. Indeed, when x G Vk, 
owing to the singularity of the kernel GK,{x,y) ai y = x, the integrand is unbounded within 
the integration domain and direct use of a standard quadratures yield inaccurate approximations. 
Thus, specialized quadrature rules must be developed and used to deal with such singular integrals. 
The case when x ^ Vk, in contrast, does not present this challenge. To effectively present these 
two contrasting scenarios, we further refine the index sets Ib and Ij by introducing target point 
dependent index sets M.b{x) = {A: G Xb \ x G Vk} and M.i{x) = {/c G X/ | a; G Vk} to rewrite 0 


as 


lC[u]{x) = 


E + E 

.fceA4s(a;) 


( 10 ) 



G^,{x, y)m{y)u{y)ujkiy) dy 


+ ' E + E )/ G^{x, y)m{y)u{y)ujk{y) dy. 

Clearly, each of the integrals above can be rewritten in the parametric coordinates 

j Gn{x,y)m{y)u{y)ujk{y)dy = jj G^{x,^k{ti,t2))(pk[u\iti,t2)Ckiti,t2) dtidt2, 

Vk [0,1]2 

where ^k[u\{ti,t 2 ) = rn{^kiti,t 2 ))u{^k{ti,t 2 ))i 0 k{$k{ti,t 2 )), and denotes the Jacobian of the 
transformation ^k- 


Now, for the cases k 0 XA.b{x) and k 0 Aii{x) in (10), G^ix, y) remains non-singular throughout 
the region of integration. Though adopting a single high-order approximation scheme for both 
scenarios is possible, we actually utilize two different quadratures to take advantage of a more 
favorable behavior of the integrands in the later case when k 0 Aij{x) where (/?fc[u] vanish to 
high order at the boundary of the integration domain and have smooth periodic extensions to 
As is well known, the trapezoidal rule exhibits super-algebraic convergence for smooth and 
periodic integrands, which we indeed employ to obtain accurate approximations in this case. In 
contrast, when k 0 M.b{x), a straightforward use of trapezoidal rule does not produce high- 
order accuracy as the integrands do not vanish at t 2 = 0. We overcome this minor difficulty 
by utilizing the trapezoidal rule for the integration with respect to ti-variable while employing 
a composite Newton-Cotes quadrature in transverse variable t 2 to achieve approximations whose 
rate of convergence directly depends on the order of Newton-Cotes used and could be enhanced 
arbitrarily as long as the smoothness of m within 11 allows it. 

As mentioned above, integrands in (10) corresponding to the cases when k G AJs(a:) and 


k G Ai[{x) are singular. To achieve rapidly convergent approximations, we rely on analytic resolu¬ 
tion of singularities through suitable changes of parametric variables and application of high-order 
quadratures to resulting smooth integrands. In addition, we also localize the region where such co¬ 
ordinate transformations are affected to a small neighborhood of the singular point using a suitable 




6 


AKASH ANAND, AMBUJ PANDEY, B. V. RATHISH KUMAR, JAGABANDHU PAUL 


smooth and compactly supported cut-off function. Indeed, we see that 

( 11 ) 

G^,{x, ^fc(t)) ■■■ dt = jj G^{x, ^fc(t)) • • • r]{t; ^^^{x)) dt+ j ^ G^{x, ik{t)) ■ ■ ■ (!-??(*; ^k^{x))) dt, 
[ 0 , 1]2 [ 0 , 1]2 [ 0 , 1]2 

where rj{-;to) is a G°^ function that it is compactly supported in a neighborhood of the point 
to £ [0,1]^ while ry = 1 in a smaller neighborhood of to- This localization of singularity as seen in 
the first integral on the right hand side of 0 brings in a two-fold benefit, namely, (a) it limits 
the relatively expensive treatment of singularity which, among other computational challenges, also 
demands an interpolation of grid-data to off-grid quadrature points to a small integration domain, 
and (b) it allows for additional speed-up in the computation of non-singular second term on the 
right hand side of ©• Note that this regular integral, of course, can be integrated to high-order 
using the same numerical quadratures that we apply in the case of a; 0 "Pfe, as explained above. 

It is straightforward to see that the direct application of the integration scheme that we described 
above leads to a computational complexity of O(iV^), where N = |T| is the size of the set T, the 
total number of quadrature points in the Nystrom scheme. However, as we outline below, one can 
improve the computational complexity of this methodology to 0{N log N) by breaking the overall 
computation of these integrals into two parts - a relatively small but specialized calculation in the 
neighborhood of singularity and remaining non-singular contributions arising from the voluminous 
bulk. Toward this, we begin by bounding the inhomogeneity H by a square cell C of side length A. 
This cell is further partitioned into identical cells Cij {i,j = 1, • • • , L) of side length H = A/L 
such that the bounding cell C contains L of these smaller cells along its sides (for an example, see 
Figure]^. We assume, without loss of generality, that for the chosen side length A, cells Cij do not 
admit inner acoustical resonance. This, of course, can be easily ensured by necessary adjustment 
in the choice of A so that —k? is not a Dirichlet eigenvalue of the Laplace operator in the cell C, 


These cells are used to break the computation of integrals appearing in equation (10) into adjacent 


ly 


and non-adjacent calculations which, as pointed out above, is primarily motivated by our desire 
for the method to have a more favorable computational complexity than the quadratic cost in the 
number of unknowns. 

We say a source point y is adjacent to a fixed point target point x G Cij , if it belongs to the set 
Af{x) defined by 

(12) Af{x) = {y \ y G Gki for some k, I satisfying |fe — i| < 1, |/ — jj < 1} . 

Obviously, we say y is non-adjacent to x ii y ^ J\f{x). Based on this, we separate the integrals 


over Vk in (10) into integrals over Vk riJ\f{x) accounting for the adjacent contributions and the 


computationally large non-adjacent contributions for which a variant of the two face FFT based 
acceleration strategy introduced in |29j is employed for efficient calculations. 


We, thus, split integrals in (10) to rewrite the expression for /C[u](a;) as 


^N(®) = ^(0 + 0 / +Y. . } Gk,{x, y)m{y)u{y)uJk{y) dy 

.keMsix) k^Mslx)) 


(13) + / } G^{x,y)m{y)u{y)u}k{y)dy, 

yk&Mi{x) k^Miix)} '’'Pkr\N{x) k&Xi'’'^kr\N{xY 
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where M{xY denotes the compliment of M{x), that is, the set T \ M{x). 

We note here that the computational cost of carrying out the adjacent calculations for all target 
points £c G T, in view of the fact that on an average each cell contain O discretization 

points, is X O (A^/L^) = O (A^^/L^). 

The cost involved in computing the contributions to inte grals coming from Vk H N{xY for all 
target points a; G T, as we explain later in Section (3.2) is given by log(LA^^/^)) + 

/LY) + /L). Thus, we get 


0(LN2 log (LN 2 



as the total cost of computing /C[u] at all Nystrom nodes. This, of course, suggests that by choosing 
the parameter L = 0(A^^/^), the computational complexity of the algorithm reduces to the desired 
0{N log N). 

3.1. Singular Integration. Recall that when the target point x belongs to the integration patch 
Vk, we break the integral as follows: 

j Gn{x,y)m{y)u{y)u}k = jj GKix,$kit))(pk[u\it)^kit)iy)r]itl^Y^ix)) dt+ 

Vk [0,1]2 

G«(a;,^fc(t))V9fc[u](t)^fc(t)(?/)(l - y{Yi-^{x)))dtV 

( 14 ) j G^{x,y)m{y)uiy)u;k{y)il-v{^f:\y)'^^k^i^)))dy■ 

Vkn^f(x)^ 


As mentioned above, we rely on the FFT-accelerator for computation of the third term in (14) 


This, however, entails that we compute this integral as a convolution which, in turn, requires that 
the integrand remains independent of localization cut-off function y. This, of course, can easily be 
achieved by ensuring that the support of is contained in M{x). The second term in 

this expression, which we refer to as local correction to the accelerated numerics, is computationally 


small. We delay a more detailed discussion on these aspects to Section (3.2). We detail the 
approximation of the singular integration 

(15) fCY'^^iuYx) = jj G^{x,$kit))Tk[u]{t)Ykit){y)r]{t]^j^\x))dt, 

[0,1]2 n supp?; 


the first term on the right hand side of equation (14), next. 


3.1.1. Singular integration over interior patches. In this section, we describe high-order quadra¬ 
ture rule for evaluation of /C^*"'^[u](®) when k G M.i{x). In this case, following [29], we use the 
localization function r]{t; (x)) = x{\t — ^(®)l/^) with an appropriate real number r > 0 and 

a smooth function y such that x = 1 in a neighborhood of 0 and x(s) = 0 for all s > 1. We extend 
the domain of integration in (15) to the disc of radius r around = (ti,^®) each target 


point X £VkF\^. Obviously, because of the fact that vanishes to high order in all directions 
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on the boundary of the integration domain, this process does not compromise on the smoothness of 
the integrand. To perform the integration, we then change to polar coordinates centered at (t*, ): 

ti = ti + pcosO, t2 = t2 + psinO. 

If we let (pk [u] {tf + p cos 0, -\- p sin 9) = 

(^kN(ti + /9cos0,tf + /?sin0)^^(t* + pcos 6 ,t 2 + psm9)p{{tf + /?cos0,tf + /3sin0); (ti 


then the integral (15) takes the form 


2n 


(16) [u\{x) = - d9 \p\GK{x,^kiti+pcos9,t2+psm9))ipk[u\{ti+pcos9,t2+psin9) dp. 


The appearance of additional factor \p\ (the Jacobian of the polar change of variables) in the 
integrand of integral (16) cancels the kernel singularity as \p\Gi.^{x,y{ti + pcos0,tf + psin9)), 


clearly, is a smooth function of p. Additionally, the cut-off function rj vanishes to high order at the 
boundary of integration interval in p variable, and therefore, a use of Trapezoidal rule yields super- 
algebraic convergence for the p-integration. Moreover, the 0-integral can also be approximated 
to high order by employing Trapezoidal rule quadrature as the corresponding integrand varies 
smoothly and periodically with respect to 9. 

While the change to polar coordinates provides a way to resolve the singularity analytically, the 
proposed application of Trapezoidal rule demands that we provide the values (pfc[u] at points outside 
of the computational grid T. This necessitates employing an efficient and accurate interpolation 
strategy for evaluation of <Pfc[u] at these newly transformed grid points. Toward this, we adapt the 
Fourier refined polynomial interpolation introduced in [29], with a suitable choice of polynomial 
degree, to retain high order accuracy while maintaining computational efficiency. 


3.1.2. Singular integration over boundary patches. When target point x lies in one of the boundary 
integration patches, say Vk with k G Aisix), we adopt a different approximation strategy to 
evaluate /C^*”'^[u](a;). The first divergence comes in the form of the choice of the cut-off function 
rj where a circular support used for interior patches becomes unsuitable for use near the physical 
boundary of the domain, that is, near t 2 = 0. In this case, therefore, following the ideas introduce 
in m, we take the localization function to have a rectangular support around ^(a;) = (tf jtf); 
given by 

V{{h,t2y, {tf,tf)) = x(|ti - tf\/r)x{\t2 - if l/r). 

We represent /C^*"^[u](£c) as 


(17) 


where 


(18) 


= j Jk{t2;x)x{\t2 - t2\/r)dt2 

0 


tf+r 


Jk{t-,x)= / Gt^{x,$,k{ti,t))pk[u\{ti,t)x{\ti-kl\/r)dti. 


tf-r 
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The evaluation of x) poses difficulties owing to singularity present in the integrand at t 
and the near singularity in the vicinity of . In order to circumvent these difficulties, we use a 
change of variable ti = tf + q{t), where the smooth invertible odd function ^j(t) satisfies 


(19) 


d'^gir) 


dr^' 


T=0 


= 0 for m = 0,M, 


(for example, £»(t) = for even, non-negative integer M). 

variable in (18) to obtain 


We now change the integration 


‘(r) 


( 20 ) 


(T ^) 


G{x,^k{tf + giT),t))ipk[u\{tf + Q{T),t)x{Q{T)/r)Q'{T) dr. 


-e bO 


This change of variable renders the integrand M times differentiable in r that are also uniformly 


bounded. In addition, in view of the factor x(^>(T)/r), the integrand in equation (20) vanishes on 
the boundary of the integration interval i.e. r = ±^“^(r), together with all of its derivatives. Thus 


the integrand in integral (20) is a smooth and periodic function which can be integrated to high- 


order by means of Trapezoidal rule. Again, this change in a variable produces a set of quadrature 
point in r that do not coincide with the computational grid on the integration patch Vk requiring 
an interpolation strategy for which we utilize the TFT refined polynomial interpolation |29| . 

The final step in high-order approximation of /C^*”^[u](£c) relates to the the computation of t 2 - 
integral in ( |17[ ). The main difficulty, here, is encountered in the form of a jump discontinuity in 
the t-derivative of 74(t; x) ai t = [IS]- To work around this, we split the integral in as 


( 21 ) 


‘'2 1 
= J Jk{t]x)dt + j Jk{t;x)dt, 


where both integrands, in principle, can be approximated to high-order by means of Q-point 
Newton-Cotes quadrature. This, however, presents a practical difficulty in the form of requir¬ 
ing at least Q-equidistant grid points in [0,tf], and [tf,l], that, of course, is not available when 
is close to either 0 or 1. The implementation of a Q-point quadrature, in such cases, requires 
values of 74(4 ®) a-t points other than the original grid points. The direct interpolation of 74(4 x), 
however, is neither accurate, (on account of non-smoothness of 74(4®) at t = tf), nor efficient 
(because of its dependence on tf). To avoid this expensive computation, we further split (21) as 
follows: 


fif—Li(Q—l)/i2 —(ii —1)^2 


/ Jk{t]x)dt+Y] 

^2 1*^2 H “^2 ( Q ~ 1)^2 ~^2 /" t ?+^2 ( Q ~ l )^2+'^2^2 

+ E/ Jk{t;x)dt+Y. 

l 2 = ± “^(^2 + 1)( Q ~1)^2 22 = 1 ^2 "^^2( Q ~1)^2 + ('^2~1)^2 

pt^-L-i{Q-l)h2-Iih2 rl 

(22) -h / Jk{t;x)dt+ / Jk{t;x)dt 

Jo Vt 2+^2 (Q~ 1)^2 + ^ 2^2 


Jk{t]x)dt 


Jk{t]x)dt 
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where /12 = I/-/V 2 denotes the mesh size in t 2 direction and numbers Li, L 2 , h, I 2 are obtained as 

1-tf 

{Q-l)ht\ ’ 


tf - Li{Q - l)/i2 

Tr, — 

r(l-tf)-L 2 (Q-l)h 2 l 

h2 

? ^2 

h2 


Li = 


fX 


{Q — l)/l2 


; L 2 = 


with [r] denoting the largest integer less than or equal to the real number r. 

Note that when lies on one of the parallel grid lines, that is, = jh 2 for some j, then last two 
integrals in (22) vanish. Moreover, the computation of the second and the fourth set of integrals 
therein, requires at the most 2 x {Q — 1) x {Q — 1) additional grid points, first {Q — 1) x {Q — 1) 
in the vicinity of 0 and the remaining in the vicinity of 1. The values, at these additional 

grid points are obtained by, first interpolating the smooth density at these extra grid lines 

followed by the integration in (20). 


Remark 3. It is clear that the strategy described above is valid when the target point a; G T happens 
to lie on one of the parallel grid lines. Obviously, the overall computational scheme does require us to 
eompute [u\{x) for target point x that do not coincide with any of the grid lines. In particular, 
this happens when the integration on a boundary pateh corresponds to a target point coming from 
an interior patch. To deal with such interactions between boundary and interior patehes, we first 
precompute ICk[u\{x) at all boundary discretization points x G ^k&XB^k o,nd then set up an FFT- 
refined polynomial interpolation schemes in the boundary region Ufcgi^'Pfc making possible aecurate 
evaluations of 1Ck[u\ at off-grid lines. 


3.2. Non-singular integration. As mentioned in the beginning of this section, for the fast com¬ 
putation of non-singular integrals 


(23) 


IC^^3[u]{x) = 



PkrWix)^^ 


G,,{x, y)m{y)u{y)u}kiy) dy, 


that appear in (13), we rely on an equivalent sources approximation strategy [29] to improve the 


overall computational efficiency. While the equivalent source methodology is given for the three 
dimensional Helmholtz kernel in [29|, we can readily adapt those ideas to obtain a variant for the 
two dimensional case, as discussed below. 

Indeed, we can accurately approximate /C^'^^[u](a;q) at each grid point aig G T at an O(N^) cost 
by employing a high order quadrature, say, given by 


(24) 


)Cr[u]{x,) = 


E 

yteX\J\f{xq) 


wiGn{xq, ye)m{yi)u{yi)ujkiye). 


Noticing that (24) is a discrete convolution whose contributing sources yi are located somewhat 


irregularly in H, the acceleration strategy seeks to replace them by a certain set of “equivalent 
sources” placed on a Cartesian grid, that produce an accurate approximation for the convolution, 
to facilitate the use of FFT for computing the discrete convolution at an 0{N log N) cost. 

More precisely, for each cell Cij, if IC^J^^j[u]{x) denotes the quantity 


^a,ij N (®) = E WiG,^{x,yi)m{ye)u{yi)ujk{yi), 
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Figure 2. An illustration of the locations Xij^i of equivalent sources (points dis¬ 
played in blue color) on the parallel faces of cells Cij. The image also depicts, for a 
target point x, the set M{x) that defines the adjacency. 


then we seek constants and defining the approximating quantity 


(25) 


such that 


N<^q 


Am), 


)d) dG{x, Xij^e) 


J\J'COll 


q=l 

is minimized for a fixed number of evaluation points Xq on the boundary of the set Af{x) 

corresponding to the cell Cij. The locations Xij^i of equivalent sources in (25) are equidistant points 
placed on the two opposite and parallel sides of the cell as shown in Figure It is known that one 
can achieve this approximation up to the prescribed accuracy 0{e), provided the number is 
chosen as ini 

It 21og(3/\/2) 

The unknown constants and \ are obtained as solution to the overdetermined linear system 
Acr^^ = b. It is important to note that, as the geometry is identical for each cell Cij, the QR 
factorization of the above matrix A need only be computed once and saved for repeated use. For 
numerical stability of the least-squares solver, we choose A straightforward counting 
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show that this process of equivalent source computation for cells requires /L^) + /L) 

operations in total. 

Clearly, for x G Ctj, the computation of 

L L i+1 i+1 

k=l 1=1 k=i—l l=j—l 

as an approximation to /Ca^^[u](») at all points on the grid can now be obtained by means of 
FFTs with a computational cost of O log(LiV^/^)). The values obtained in this manner 

provide accurate approximations for non-adjacent non-singular interactions at points Xij^i on the 
boundary of cells Cij. Finally, to obtain these values at a true source location, say for an £c G Cij, we 
solve the free space Helmholtz equation within Cij, with Dirichlet boundary boundary data coming 
from ICa‘^‘^’^‘^[u]{xij^i). To efficiently obtain solutions to these well-posed Dirichlet boundary value 
problems, we utilize a discretized plane wave expansion of the form [30] 

(26) )Cl^'i'^'i[u]{x) « '^jje^piinde.x), 

£=1 


where the unit vectors dj sample the surface of unit disc with sufficient degree of uniformity. Use of 
this approach is motivated by the spectral convergence of the above expansion with respect to the 
number of unit vectors di used and the fact that the wave expansion coefficients 'ji can be obtained 
as solutions to overdetermined linear systems of the form By = (3ij where the matrix B remains 
unchanged for each cell, again, owing to the identical geometry of cells CijS. Thus, ICa'^^’^‘^[u\{x) can 
be evaluated at all true source locations x G Cij at a computational cost of /L^)+0{N^^'^/L) 

and the overall cost of evaluation of nonadjacent non-singular interactions, therefore, stands at 
O (LiVV2 log(LAfV2)) + 0{N^/^/L^) + 0{N^/‘^/L). 

3.3. Numerical Results. In this section, we present a series of numerical examples to exemplify 
the effectiveness of the two dimensional algorithm discussed in previous subsections in terms of 
its numerical accuracy and computational efficiency. Be begin by noting that, in all convergence 
tables given below, we use x ni x n 2 -|- O'! x mi x m 2 to announce the underlying computational 
grid implying that a pi number of boundary patches, each with ni x n 2 discretization points, and 
a q 2 number of interior patches, each with mi x m 2 points, constituted the Nystrom grid T. The 
relative error (in the near field) reported in tables are obtained as 


£00 


max 

l<i<7V 


^exact 


(®i) -U^PP™^(£Ci)| 


max 

l<i<N 


£2 



exact 


i=l 


{Xi) - u'^PP™^(®i 


2\ 


1 

2 


N 
2 = 1 


We write “Order” to denote the numerical order of convergence of the approximation. All conver¬ 
gence studies correspond to the scattering of incident plane waves propagating along the x-axis. 
As our approximation scheme for the integral operator is based on the use of spectrally accurate 
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Grid Size 

Unknown 

Iteration 

Number 

L'^ 

L°° 

£2 

Order 

^OO 

Order 

2x3x9-|-lx9x9 

135 

10 

7.04e-01 

- 

6.21e-01 

- 

2 X 5 X 17-M X 17x17 

459 

11 

2.68e-01 

1.39e-h00 

2.73e-01 

1.18e-F00 

2x9x33-Mx 33x33 

1683 

12 

3.92e-02 

2.78e-h00 

4.77e-02 

2.52e-F00 

2 X 17x 65-bl X 65x65 

6435 

15 

3.38e-03 

3.54e-h00 

4.14e-03 

3.52e-F00 

2 X 33 X 129 -h 1 X 129x129 

25155 

18 

1.47e-04 

4.52e-h00 

2.37e-04 

4.13e-F00 


Table 1 . Convergence study: plane wave scattering by a disc of acoustic size na = 
10 and refractive index n = \/2 when 3-point Newton-Cotes quadrature employed 
for t 2 -integration over boundary patches. 


Grid Size 

Unknown 

Iteration 

Number 


LOO 

£2 

Order 

^OO 

Order 

2x5x9-|-lx9x9 

171 

10 

7.09e-01 

- 

5.96e-01 

- 

2 X 9 X 17-M X 17 X 17 

595 

11 

1.22e-01 

2.54e-k00 

1.27e-01 

2.24e-F00 

2 X 17 X 33 -h 1 X 33 X 33 

2211 

15 

8.78e-03 

3.80e-k00 

1.03e-02 

3.62e-F00 

2 X 33 X 65 -M X 65 X 65 

8515 

18 

2.84e-04 

4.95e-h00 

3.35e-04 

4.94e-F00 

2 X 65 X 129 -h 1 X 129 X 129 

33411 

23 

7.76e-06 

5.19e-h00 

8.53e-06 

5.30e-F00 


Table 2. Convergence study: plane wave scattering by a disc of acoustic size Ka = 
10 and refractive index n = \/2 when 5-point Newton-Cotes quadrature employed 
for t 2 -integration over boundary patches. 


Trapezoidal rule and Q-point Newton-Cotes quadrature, the expected convergence rates are essen¬ 
tially that of Newton-Cotes, that is, either Q or Q + 1 depending on whether Q is even or odd, 
provided other parameters, such as the degree of the polynomial in the interpolators, the number 
M for the change of variable g, etc. are chosen favorably. We demonstrate the increase in the 
convergence rate with respect to increasing Q through two different sets of experiments, one using 
a 3-points while other using 5-points Newton-Cotes quadrature. 

Example 3.1. {Convergence study for a simple scatterer) 

As a first exercise, we compute the scattering by a penetrable disc for which a comparison with 
exact solution is possible. For this scattering configuration, we present a convergence study, in 
Table when the acoustic size {na, a being the diameter of the inhomogeneity) of the scatterer 
is 10 and its refractive index is n = \/2. The numerical solutions in this study are computed for 
several discretization levels while employing a 3-point Newton-Cotes quadrature for t 2 -integration 
in boundary patches. In this case, we clearly see the methodology achieve the expected rate of 
convergence. To show an enhanced order of convergence, we repeat the same experiment but use a 
5-point Newton-Cotes quadrature in the computation of the integral operator. The corresponding 
results are presented in Table where, as expected, we see an increased rate of convergence. 

We display, in Figure the numerical solution obtained using our methodology corresponding 
to an experiment with na = 20 and n = \/2 where the incidence wave impinges on the obstacle at 
45 degree angle with positive x-axis. This computation required 57 GMRES iterations to produce 
an accuracy of 0.006% in the max norm. 
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(a) Disc 



(b) Plane wave incidence, 5R(u*) 




(d) Real part of the total field, 5R(u) 


Figure 3. Scattering of a plane wave exp(i/td • r) with k = 10, d = (l/\/2, l/\/2) 
by a penetrable unit disc. A computational grid of size 2 x 65 x 129 + 1 x 129 x 129 
and 5-point Newton-Cotes quadrature is used to obtain an accuracy of 0.006%. 


Example 3.2. (Computational Efficiency) We next demonstrate a tempered growth in computa¬ 
tional cost of the proposed accelerated algorithm while maintaining a fixed computational error by 
comparing our approximate integral operator against the continuous operator (§. Toward this, 
with every doubling of the computational grid size, we double the wave number k to keep the num¬ 
ber of points per wavelength unchanged, thus, fixing the accuracy level. For this set of experiments, 
we again use a disc scatterer with a constant refractive index re = 2. The results in Table show 
that, while for small values of N, the cost of non-accelerated algorithm is comparable to that of the 
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Ka 

Grid Size 

£2 

*^00 

Time(s)/Iteration 

accelerated 

Non-Accelerated 

4 

2 X 16 X 32 -Fl X 32 X 32 

1.18e-03 

3.15e-03 

2.57 

3.93 

8 

2 X 32 X 64 -|- 1 X 64 X 64 

3.90e-04 

5.67e-04 

8.53 

39.56 

16 

2 X 64 X 128 -h 1 X 128 X 128 

7.90e-04 

7.85e-04 

32.8 

626.47 

32 

2 X 128 X 256 -h 1 X 256 x 256 

8.64e-05 

9.05e-05 

122.31 

10785.9 


Table 3. Computational cost for accelerated and non-accelerated algorithm. 


kH 

Equivalent Source/cell 

£2 

^OO 

2 

4x4 

5.10e-05 

1.30e-04 

4 

8x4 

3.05e-09 

1.94e-08 

8 

13 X 4 

3.00e-12 

6.70e-12 

16 

30x 4 

1.74e-12 

1.14e-12 


Table 4. Accuracy of acceleration when wavenumber k increases but number of 


point per wavelength is hxed. 


Equivalent Source/cell 

£2 

^OO 

4x4 

1.60e-02 

4.66e-02 

6x4 

4.14e-05 

1.68e-04 

8x4 

3.53e-08 

1.85e-07 

10 X 4 

1.60e-10 

1.48e-09 

12 X 4 

2.86e-12 

2.82e-ll 


Table 5. Accuracy of acceleration for hxed wavenumber k = Svr 


accelerated version, as problem size increase, there is a substantial gain in terms of computational 
cost. For instance, for the scatterer of size no = 32, accelerated computations are 88 times faster 
than its non-accelerated counterpart. In particular, the time for accelerated computation in Table 
[^exhibit the growth according to the computational complexity of 0(iVlog A^). 

Example 3.3. {Accuracy of non-singular non-adjacent interactions) 

The previous examples demonstrate the high-order convergence and efficient computational times 
for the proposed methodology. We, in particular, noted signihcant gains for the accelerated scheme 
in terms of computational cost. On the other hand, though Table does show that the accuracy 
levels of the numerical method are maintained even when frequency of wave oscillations increase, 
provided we ensure a hxed number of discretization points per wavelength, an interesting picture 
emerges when we look approximations coming from the accelerator alone. To see this, we compare 
the values of [u] against that of [^] increasing levels of frequency while the number of 

equivalent sources remained unchanged. As seen in Table approximations converge rapidly with 
increasing wavenumber. This observation, in turn, is utilized in a more effective load balancing 
between adjacent and non-adjacent calculations to achieve favorable computational cost. 
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We also include, in Table a similar study, this time keeping the wavenumber constant (k = Svr) 
while increasing the number of equivalent sources, where we again see rapid improvement in 
accuracy levels. Tables and clearly demonstrate the spectral accuracy of acceleration strategy 
for small as well as large cell sizes H. Indeed, error introduced as a result of acceleration procedure 
are typically much smaller in comparison with other sources of error within the overall algorithm, 
and have no impact on convergence rates, as confirmed by the numerical studies in this section. 



(a) Bean 



(c) Real Part of the scattered field, IR(id) 



(b) The incident plane wave, 5R(u*) 



(d) The real part of the total field, 5R(u) 


Figure 4. Scattering by penetrable bean shape scatterer, At frequency k = W,d = 
(1, 0), Grid Size 2 x 65 x 129 + 1 x 129 x 129, Newton-Cotes hve point quadrature 
is used for transverse integration over boundary patches. 
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Grid Size 

Unknown 

Iteration 

Number 


L 

30 

£2 

Order 

^00 

Order 

2x3x94-1x9x9 

135 

4 

1.95e-400 

- 

1.82e-400 

- 

2x5xl7-4lx 17x17 

459 

11 

3.89e-01 

2.33e-400 

3.84e-01 

2.25e-400 

2x9x33-4lx 33x33 

1683 

15 

3.88e-02 

3.32e-400 

4.07e-02 

3.24e-400 

2 X 17 X 65-4 1 X 65x65 

6435 

20 

4.10e-03 

3.24e-400 

4.83e-03 

3.08e-400 

2 X 33 X 129-4 1 X 129x129 

25155 

25 

2.28e-04 

4.17e-400 

4.55e-04 

3.41e-400 


Table 6. Convergence for the bean shape scatterer with k = 5 and n = y/2, 
when 3-point Newton-Cotes quadrature employed for boundary patch integration in 
transverse integration. 


Example 3.4. {Convergence study: complex scattering media) 

The methodology presented in this text can, of course, be applied to acoustic scattering calcula¬ 
tions from penetrable media with complicated geometries as well as variable material properties. 

We begin with an example that demonstrates adaptability and applicability of our algorithm 
in dealing with scatterers that have relatively complex geometrical description. Toward this, we 
consider scattering by a penetrable bean shaped scatterer, as shown in Figure]^ For these scattering 
computations, we use an incident plane wave exp{iKx) with ko = 10. The refractive index of the 
medium is taken to be n = \/2. We use a numerical solution at a fine grid for convergence study 
comparisons. The first sets of computations, for which the results are presented in Table employ 
3-point Newton-Cotes quadrature for t 2 -integration over boundary patches. As expected, we see a 
convergence rate of order 4 in the solution. We repeat the experiment with 5-points Newton-Cotes 
quadrature and report the results in Table where we again see an enhanced order of convergence. 

Of course, our algorithm is not restricted to consideration of constant or piecewise constant 
material properties. In order to demonstrate this fact, in our final experiment, we carried out 
scattering computations for a variable refractive index n given by 

n{x) = sin(7rxi) cos('7rx2) for x = (xi, X 2 ) G fl. 

A plane wave incidence with k = 2 is used and results corresponding to three and five point 
Newton-Cotes quadrature are presented in table Q and Q respectively. Convergence study of 


Grid Size 

Unknown 

Iter 


LOO 

£2 

Order 

^00 

Order 

2x5x9-41x9x9 

171 

12 

1.85e-400 


1.31e-400 

- 

2 X 9 X 17-4 1 X 17 X 17 

595 

16 

1.56e-01 

3.57e-400 

3.35e-01 

1.96e-400 

2 X 17 X 33 -4 1 X 33 X 33 

2211 

18 

1.90e-02 

3.04e-400 

2.94e-02 

3.51e-400 

2 X 33 X 65 -4 1 X 65 X 65 

8515 

24 

6.13e-04 

4.95e-400 

6.36e-04 

5.53e-400 

2 X 65 X 129-4 1 X 129 X 129 

33411 

30 

1.17e-05 

5.71e-400 

1.30e-05 

5.61e-400 


Table 7. Convergence for the bean shape scatterer with k = 5 and n = y/2, 
when 5-point Newton-Cotes quadrature employed for boundary patch integration in 
transverse integration. 
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Grid Size 

Unknown 

Iteration 

Number 


L 

30 

£2 

Order 

^OO 

Order 

2x3x9-|-lx9x9 

135 

3 

l.Ole-hOO 

- 

1.25e-h00 

- 

2x5xl7-Mx 17x17 

459 

4 

7.60e-02 

3.74e-h00 

9.11e-02 

3.77e-h00 

2x9x33-Mx 33x33 

1683 

5 

7.71e-03 

3.30e-h00 

8.36e-03 

3.45e-h00 

2 X 17 X 65-hi X 65x65 

6435 

6 

6.56e-04 

3.55e-h00 

6.36e-04 

3.72e-h00 

2 X 33 X 129-hi X 129x129 

25155 

7 

2.02e-05 

5.02e-h00 

5.79e-05 

3.46e-h00 


Table 8. Convergence for the 
sin(7rx) cos(7ry) for x = {x, y) G 
ployed for boundary patch integration in transverse integration. 


bean shape scatterer with k = 2 and n = 
n, when 3-point Newton-Cotes quadrature em- 


Grid Size 

Unknown 

Iteration 

Number 


LOO 

£2 

Order 

^OO 

Order 

2x5x9-hlx9x9 

171 

3 

5.25e-02 


1.69e-01 

- 

2 X 9 X 17 -h 1 X 17 X 17 

595 

4 

9.12e-03 

2.53e-h00 

3.34e-02 

2.34e-h00 

2 X 17 X 33 -h 1 X 33 X 33 

2211 

5 

7.99e-04 

3.51e-h00 

3.14e-03 

3.41e-h00 

2 X 33 X 65 -h 1 X 65 X 65 

8515 

7 

4.51e-05 

4.15e-h00 

8.14e-05 

5.27e-h00 

2 X 65 X 129 -h 1 X 129 X 129 

33411 

8 

6.14e-07 

6.20e-h00 

1.23e-06 

6.04e-h00 


Table 9. Convergence for the bean shape scatterer with k = 2 and n = 
sin(7rx) cos(7ry) for x = {x, y) G O, when 5-point Newton-Cotes quadrature em¬ 
ployed for boundary patch integration in transverse integration. 


these tables clearly support that the high-order nature of our algorithm, as expected, continues to 
hold for scattering media with variable material properties. 

For the sake of pictorial visualization, results of scattering computation for ko = 20, n = \/2 are 
shown in Figure Q for a plane wave incidence traveling along positive x-axis. The plots visualize 
a solution obtained at the end of 82 iteration of GMRES when residual had reached 10“^. 


4. An extension to three dimensions 

Though this text is primarily dedicated to presenting a fast and accurate computational strategy 
to solve the inhomogeneous acoustic scattering problem in two dimensions, as we mentioned in 
Section 2, this methodology has a straightforward extension that allows for numerical solution 
of the corresponding three dimensional counterpart. The difficulty in high-order evaluation of 
the integral operator in three dimensions are largely analogous to those that appear in the two 
dimensional setting. Given that the main algorithmic steps remain unchanged, in this section, 
we avoid much of the repetitions and only briefly highlight some of the salient points underlying 
the extension. We then present a numerical verification of the fact that this methodology, indeed, 
produce rapidly converging numerical solutions to the inhomogeneous scattering problems. 

We begin by recalling that, as in the integral /C[u](a;) is written as a sum of integrals over 
interior and boundary patches. In three dimensions, each of these integrals take the following form 
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Grid Size 

Unknown 

Iter 


L°° 

£2 

Order 

^00 

Order 

2x5x5x5-|-lx5x5x5 

375 

2 

2.03e-01 

- 

2.66e-01 

- 

2x9x9x9-hlx9x9x9 

2187 

2 

5.21e-02 

1.96e-F00 

2.48e-01 

l.Ole-01 

2x17x17x17-^1x17x17x17 

14739 

4 

5.02e-03 

3.38e-F00 

3.28e-02 

2.92e-h00 

2 X 33 X 33 X 33 -|- 1 X 33 X 33 X 33 

107811 

6 

1.65e-04 

4.93e-F00 

7.59e-04 

5.44e-h00 

2 X 65 X 65 X 65 -h 1 X 65 X 65 X 65 

823875 

8 

4.55e-06 

5.18e-h00 

1.61e-05 

5.55e-h00 


Table 10. Convergence study for scattering of plane wave by a penetrable spherical 
shape scatterer with wave number k = 2 and refractive index n = \/2, when 5-point 
Newton-Cotes quadrature employed for boundary patch integration in transverse 
integration. 


when expressed in the parametric space variables: 


Gi^{x,y)m{y)u{y)ujk{y) dy = // / dtidt 2 dt 3 . 


Vk 


[0,1]3 


As in the two dimensions, for the cases k 0 A4b{x) and k 0 A4j{x), the kernel GK,{x,y) remains 
non-singular within the integration region and, therefore, corresponding integrals can be approx¬ 
imated using high-order quadratures. Again, when k 0 the integrands have smooth and 

periodic extension to and the trapezoidal rule yields approximations with super-algebraic con¬ 
vergence. For the case when k 0 Mb{x), on the other hand, high order can be attained by simply 
employing spectrally accurate Trapezoidal rule for planar integrations (with respect to ^ 1 ,^ 2 ) and 
a high-order composite Newton-Cotes quadrature for the integration in the transverse variable ts. 

For the cases k G M.b{x) and k G M.j{x), as before, high-order evaluation of the integrals 
becomes difficult owing to the kernel singularity ak, x = y. Toward dealing with this in bound¬ 
ary patches, following [T3], we change to polar coordinates (p, 0) centered around 
the projection of on to the ts-integration plane, which, upon employing an accompanying 

polynomial change of variable in p, provides an effective resolution of difficulties arising out of the 
kernel singularity. Subsequent application of Trapezoidal rule in ti, t 2 variables and a composite 
Newton-Cotes quadrature for t 3 integration yields accurate results. For the singular integral over 
interior patches, on the other hand, we change to spherical coordinates (p, 6, cj)) around target point 
This, in turn, yields smooth integrands, where accurate integrations can be affected by 
employing Trapezoidal rule in p, cj) variables and Clenshaw-Curtis quadrature in 0 variable. 

The integration scheme can again be accelerated by a suitable use of two face equivalent source 
approximations on Cartesian grids. Just as in two dimensions, this strategy of employing three 
dimensional FFTs for approximating the convolution, further reduces the evaluation cost of the 
non-singular non-adjacent interactions. In this case, the algorithm exhibits the theoretical compu¬ 
tational complexity of 0{N\ogN) with respect to the grid size N, provided we choose L = 
when number of sub-cubes are used to cover the inhomogeneity. 

In order to demonstrate the high-order convergence of the three dimensional scheme, we present 
results from a plane wave scattering by a penetrable spherical inhomogeneity with constant refrac¬ 
tive index, n = y/2. The convergence study, given in Table 10, that corresponds to the incident 
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(a) Spherical shape scatterer (b) Intensity of the total field (|u|^) 


Figure 5. The figure of the left depicts a the discretization on the outer surface 
of the scatterer coming from a typical volumetric computational grid. The figure of 
the right presents a visualization of the total field intensity. 


plane wave u® = exp(zK 2 ;) with wavenumber k = 2, clearly exhibits rapid convergence of numerical 
solutions. 

As a final example, we visualize, in Figurej^ the results of a scattering computation corresponding 
to a spherical obstacle of acoustical size na = 10 and refractive index n = 2. For this simulations, a 
relative error of order 10“^ is achieved at a computational grid of size 2 x 33 x 65 x 65+1 x 65 x 65 x 65. 


5. Conclusions 

In this paper, we discuss a fast high-order method for scattering of acoustic waves by penetrable 
inhomogeneous media which require 0{NlogN) computational cost for each iteration of iterative 
linear system solvers. Rapidly convergent approximations are obtained through a use of specialized 
quadrature rule while a reduced computational cost is achieved by utilizing a strategy based on 
equivalent sources approximation on the Cartesian grids. We present a series of numerical exper¬ 
iments to demonstrate its performance in terms of computational efficiency as well as numerical 
accuracy. We emphasize that this algorithm is designed to retain high-order convergence even 
when material properties jump across the scattering interface and can be employed for simulating 
scattering of acoustic waves by geometrically complex structures. For instance, in our numerical 
experiments where function m = 1 — is discontinuous across the material interface which, of 
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course, limits the global smoothness of the total field u, (in fact, u G C 1(M3) [31]), we still ob¬ 
tain rapidly convergent approximations. This, of course, becomes possible as a result of carefully 

avoiding integrating across the material interface. 
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